Log-Laplace Distribution — a log-symmetric, heavy-tailed model on \((0,\infty)\)#
The log-Laplace distribution (SciPy: scipy.stats.loglaplace) is a positive distribution whose logarithm is Laplace (double-exponential).
It’s a useful choice when you want a distribution on positive magnitudes with:
symmetry on a multiplicative (log) scale,
power-law tails (heavier than lognormal), and
a simple inverse-CDF sampler.
What you’ll learn#
How to define the log-Laplace PDF/CDF and connect it to the Laplace distribution via a change of variables.
Which moments exist (and why some don’t), including closed forms for mean/variance/skewness/kurtosis when they do.
A NumPy-only inverse-CDF sampler and how to validate it.
Simple inference for the shape parameter: MLE, exact confidence intervals, and a conjugate Bayesian update.
import platform
import numpy as np
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots
import scipy
from scipy import stats
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook") # CKC convention
np.set_printoptions(precision=5, suppress=True)
rng = np.random.default_rng(7)
print("Python", platform.python_version())
print("NumPy", np.__version__)
print("SciPy", scipy.__version__)
Python 3.12.9
NumPy 1.26.2
SciPy 1.15.0
1) Title & Classification#
Name: Log-Laplace distribution (
loglaplace)Type: Continuous
Support (standard): \(x \in (0,\infty)\)
Parameter space (standard): shape parameter \(c>0\)
SciPy uses an additional location/scale transform:
where scale > 0 and the support becomes \(x > \text{loc}\).
We write (standard form):
2) Intuition & Motivation#
2.1 What it models#
A log-Laplace random variable can be written as
So additive Laplace noise in log-space becomes multiplicative, heavy-tailed noise on the original scale.
A key symmetry is:
In fact, \(X\) and \(1/X\) have the same distribution.
2.2 Typical real-world use cases#
Positive data with occasional extreme values: file sizes, response times, claim sizes.
Multiplicative error models: \(Y = \theta \cdot X\) where \(\log X\) has Laplace-like “spiky + outlier-prone” behavior.
Robust alternatives in log-space: if \(\log Y\) is better modeled by Laplace than Normal (heavier tails), then \(Y\) may be log-Laplace rather than lognormal.
2.3 Relations to other distributions#
Laplace: \(\log X \sim \mathrm{Laplace}(0, 1/c)\).
Lognormal: \(\log X\) Normal vs Laplace (log-Laplace has sharper peak and heavier tails).
Pareto-like tails: for \(x\to\infty\), \(f(x)\propto x^{-(c+1)}\).
Invariance: \(X \stackrel{d}{=} 1/X\) (log-symmetry).
3) Formal Definition#
3.1 PDF#
The standard log-Laplace density can be written in a compact way:
Equivalently, as a piecewise power law (useful for intuition):
3.2 CDF#
3.3 Quantile function (inverse CDF)#
This is especially convenient for sampling:
def _validate_c(c: float) -> float:
c = float(c)
if c <= 0:
raise ValueError("c must be > 0")
return c
def loglaplace_pdf_std(x: np.ndarray, c: float) -> np.ndarray:
"""Standard log-Laplace PDF on (0, inf) with shape c>0."""
x = np.asarray(x, dtype=float)
c = _validate_c(c)
out = np.zeros_like(x, dtype=float)
mask = x > 0
xm = x[mask]
# Stable expression: f(x)= (c/(2x))*exp(-c*|log x|)
out[mask] = 0.5 * c / xm * np.exp(-c * np.abs(np.log(xm)))
return out
def loglaplace_logpdf_std(x: np.ndarray, c: float) -> np.ndarray:
"""Standard log-Laplace log-PDF."""
x = np.asarray(x, dtype=float)
c = _validate_c(c)
out = np.full_like(x, -np.inf, dtype=float)
mask = x > 0
xm = x[mask]
out[mask] = np.log(c) - np.log(2.0) - np.log(xm) - c * np.abs(np.log(xm))
return out
def loglaplace_cdf_std(x: np.ndarray, c: float) -> np.ndarray:
"""Standard log-Laplace CDF."""
x = np.asarray(x, dtype=float)
c = _validate_c(c)
out = np.zeros_like(x, dtype=float)
mask = x > 0
xm = x[mask]
out[mask] = np.where(
xm < 1.0,
0.5 * np.power(xm, c),
1.0 - 0.5 * np.power(xm, -c),
)
return out
def loglaplace_ppf_std(p: np.ndarray, c: float) -> np.ndarray:
"""Standard log-Laplace inverse CDF."""
p = np.asarray(p, dtype=float)
c = _validate_c(c)
if np.any((p <= 0) | (p >= 1)):
raise ValueError("p must be in (0, 1)")
return np.where(
p < 0.5,
np.power(2.0 * p, 1.0 / c),
np.power(1.0 / (2.0 * (1.0 - p)), 1.0 / c),
)
def loglaplace_rvs_numpy(
c: float, size: int | tuple[int, ...], rng: np.random.Generator
) -> np.ndarray:
"""NumPy-only sampling via inverse CDF (standard form)."""
c = _validate_c(c)
u = rng.random(size)
# Keep u in (0,1) to avoid returning exactly 0 or inf from floating endpoints.
u = np.clip(u, np.finfo(float).tiny, 1.0 - np.finfo(float).eps)
return np.where(
u < 0.5,
np.power(2.0 * u, 1.0 / c),
np.power(1.0 / (2.0 * (1.0 - u)), 1.0 / c),
)
def loglaplace_pdf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
"""Location/scale version: X = loc + scale*Y, Y ~ LogLaplace(c)."""
x = np.asarray(x, dtype=float)
c = _validate_c(c)
scale = float(scale)
if scale <= 0:
raise ValueError("scale must be > 0")
z = (x - loc) / scale
return loglaplace_pdf_std(z, c) / scale
def loglaplace_logpdf(
x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0
) -> np.ndarray:
x = np.asarray(x, dtype=float)
c = _validate_c(c)
scale = float(scale)
if scale <= 0:
raise ValueError("scale must be > 0")
z = (x - loc) / scale
return loglaplace_logpdf_std(z, c) - np.log(scale)
def loglaplace_cdf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
x = np.asarray(x, dtype=float)
c = _validate_c(c)
scale = float(scale)
if scale <= 0:
raise ValueError("scale must be > 0")
z = (x - loc) / scale
return loglaplace_cdf_std(z, c)
def loglaplace_ppf(p: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
c = _validate_c(c)
scale = float(scale)
if scale <= 0:
raise ValueError("scale must be > 0")
return loc + scale * loglaplace_ppf_std(p, c)
# Quick cross-check against SciPy (standard form: loc=0, scale=1)
x_test = np.array([1e-3, 0.2, 0.9, 1.0, 2.0, 10.0])
c_test = 2.5
dist = stats.loglaplace(c_test)
print("max |pdf - scipy|:", np.max(np.abs(loglaplace_pdf_std(x_test, c_test) - dist.pdf(x_test))))
print("max |cdf - scipy|:", np.max(np.abs(loglaplace_cdf_std(x_test, c_test) - dist.cdf(x_test))))
max |pdf - scipy|: 3.7947076036992655e-19
max |cdf - scipy|: 0.0
4) Moments & Properties#
A convenient fact is that the log-Laplace has power-law tails, so not all moments exist.
4.1 Raw moments#
For real \(k\) with \(-c < k < c\),
In particular:
Mean exists iff \(c>1\): $\(\mathbb{E}[X] = \frac{c^2}{c^2 - 1}.\)$
Second moment exists iff \(c>2\): $\(\mathbb{E}[X^2] = \frac{c^2}{c^2 - 4}.\)$
4.2 Variance#
For \(c>2\):
4.3 Skewness and kurtosis#
These require higher moments:
Skewness exists for \(c>3\):
Excess kurtosis exists for \(c>4\):
4.4 MGF / characteristic function#
Because the right tail behaves like \(x^{-(c+1)}\), the moment generating function
diverges for any \(t>0\) (too much mass in the far right tail for an exponential weight).
For \(t<0\), the Laplace transform exists and can be written using incomplete gamma functions:
The characteristic function \(\varphi_X(\omega)=\mathbb{E}[e^{i\omega X}]\) exists for all real \(\omega\) (bounded integrand), and admits an analogous representation with complex arguments.
4.5 Entropy#
Using \(\log X \sim \mathrm{Laplace}(0,1/c)\) and the entropy change-of-variables rule,
(Here \(h\) is differential entropy.)
4.6 Other handy properties#
Median: \(\mathrm{median}(X)=1\).
Log-moments: \(\mathbb{E}[\log X]=0\) and \(\mathrm{Var}(\log X)=2/c^2\).
Reciprocal invariance: \(X \stackrel{d}{=} 1/X\).
def loglaplace_raw_moment(k: float, c: float) -> float:
"""Return E[X^k] for the standard log-Laplace, when it exists (|k|<c)."""
c = _validate_c(c)
k = float(k)
if not (-c < k < c):
return np.inf
return (c * c) / (c * c - k * k)
def loglaplace_mean(c: float) -> float:
c = _validate_c(c)
return loglaplace_raw_moment(1.0, c) if c > 1 else np.nan
def loglaplace_variance(c: float) -> float:
c = _validate_c(c)
if c <= 2:
return np.nan
c2 = c * c
return c2 * (2 * c2 + 1) / ((c2 - 4) * (c2 - 1) ** 2)
def loglaplace_skewness(c: float) -> float:
c = _validate_c(c)
if c <= 3:
return np.nan
c2 = c * c
num = 2 * (15 * c2 * c2 + 7 * c2 + 2) * np.sqrt(c2 - 4)
den = c * (c2 - 9) * (2 * c2 + 1) ** 1.5
return num / den
def loglaplace_excess_kurtosis(c: float) -> float:
c = _validate_c(c)
if c <= 4:
return np.nan
c2 = c * c
num = 6 * (
2 * c**10 + 138 * c**8 - 615 * c**6 - 449 * c**4 - 132 * c2 - 24
)
den = c2 * (c2 - 16) * (c2 - 9) * (2 * c2 + 1) ** 2
return num / den
def loglaplace_entropy(c: float) -> float:
c = _validate_c(c)
return 1.0 + np.log(2.0 / c)
for c in [0.8, 1.5, 2.5, 5.0]:
print(
f"c={c:>4}: mean={loglaplace_mean(c)}, var={loglaplace_variance(c)}, "
f"skew={loglaplace_skewness(c)}, excess_kurt={loglaplace_excess_kurtosis(c)}"
)
c= 0.8: mean=nan, var=nan, skew=nan, excess_kurt=nan
c= 1.5: mean=1.8, var=nan, skew=nan, excess_kurt=nan
c= 2.5: mean=1.1904761904761905, var=1.3605442176870748, skew=nan, excess_kurt=nan
c= 5.0: mean=1.0416666666666667, var=0.10540674603174603, skew=3.0046141326124665, excess_kurt=40.717785467128024
5) Parameter Interpretation#
The shape parameter \(c\) is best understood via the log transform:
Larger \(c\) means smaller Laplace scale \(1/c\) in log-space, so \(X\) concentrates more tightly around \(1\).
Smaller \(c\) means heavier tails on the original scale (more extreme small and large values).
A quick rule of thumb from the tail:
So \(c\) acts like a Pareto tail index on the right tail.
6) Derivations#
6.1 Expectation (and when it exists)#
Using the piecewise PDF and splitting at 1:
The first integral is finite when \(k>-c\); the second is finite when \(k<c\). Evaluating gives:
Setting \(k=1\) shows the mean exists iff \(c>1\).
6.2 Variance#
For \(c>2\) we have \(\mathbb{E}[X]\) and \(\mathbb{E}[X^2]\), so
which simplifies to
6.3 Likelihood (iid sample) and MLE#
For data \(x_1,\dots,x_n\) with \(x_i>0\) in the standard model (only \(c\) unknown),
The log-likelihood is
Differentiating and setting to zero:
Interpretation: since \(\log X\) is Laplace, \(|\log X|\) is Exponential(rate \(c\)), so inference on \(c\) reduces to inference on an exponential rate.
def loglaplace_loglik_c(x: np.ndarray, c: float) -> float:
"""Log-likelihood for the standard model (c only)."""
x = np.asarray(x, dtype=float)
c = _validate_c(c)
if np.any(x <= 0):
return -np.inf
n = x.size
s = np.sum(np.abs(np.log(x)))
return n * np.log(c) - n * np.log(2.0) - np.sum(np.log(x)) - c * s
def loglaplace_mle_c(x: np.ndarray) -> float:
x = np.asarray(x, dtype=float)
if np.any(x <= 0):
raise ValueError("all observations must be > 0")
s = np.sum(np.abs(np.log(x)))
if s <= 0:
raise ValueError("degenerate sample: sum |log x| must be > 0")
return x.size / s
# Quick check: simulate and verify the MLE is close to truth
c0 = 2.5
x = loglaplace_rvs_numpy(c0, size=10_000, rng=rng)
print("c0 =", c0)
print("c_hat =", loglaplace_mle_c(x))
c0 = 2.5
c_hat = 2.4890218861144686
7) Sampling & Simulation#
7.1 Inverse CDF method (NumPy-only)#
Draw \(U\sim\mathrm{Unif}(0,1)\) and transform using the quantile function:
This is exactly what loglaplace_rvs_numpy implements.
7.2 Alternative viewpoint (log transform)#
Sample \(Z\sim\mathrm{Laplace}(0,1/c)\) and set \(X=\exp(Z)\). This is conceptually helpful, but the inverse-CDF above is simpler to implement without relying on a Laplace RNG.
8) Visualization#
We’ll visualize:
the PDF and CDF for different \(c\),
Monte Carlo samples (including a log-scale view).
# Grids that resolve both sides around x=1
x_left = np.logspace(-3, 0, 400, endpoint=False)
x_right = np.logspace(0, 3, 400)
xgrid = np.unique(np.concatenate([x_left, x_right]))
c_values = [0.8, 1.5, 3.0, 7.0]
fig = make_subplots(rows=1, cols=2, subplot_titles=("PDF", "CDF"))
for c in c_values:
fig.add_trace(
go.Scatter(x=xgrid, y=loglaplace_pdf_std(xgrid, c), mode="lines", name=f"c={c}"),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(
x=xgrid,
y=loglaplace_cdf_std(xgrid, c),
mode="lines",
name=f"c={c}",
showlegend=False,
),
row=1,
col=2,
)
fig.update_xaxes(title="x", type="log", row=1, col=1)
fig.update_xaxes(title="x", type="log", row=1, col=2)
fig.update_yaxes(title="density", row=1, col=1)
fig.update_yaxes(title="F(x)", row=1, col=2)
fig.update_layout(width=980, height=380, legend_title_text="shape c")
fig.show()
# Monte Carlo: samples on x-scale and log-scale
c_mc = 2.5
n_mc = 50_000
x_samp = loglaplace_rvs_numpy(c_mc, size=n_mc, rng=rng)
z_samp = np.log(x_samp)
zgrid = np.linspace(np.quantile(z_samp, 0.001), np.quantile(z_samp, 0.999), 500)
laplace_pdf = 0.5 * c_mc * np.exp(-c_mc * np.abs(zgrid)) # Laplace(0, 1/c_mc)
fig = make_subplots(
rows=1,
cols=2,
subplot_titles=("Samples on original scale (log x-axis)", "log(samples) vs Laplace(0,1/c)"),
)
# Original-scale samples (histogram) with theoretical PDF overlay
fig.add_trace(
go.Histogram(x=x_samp, nbinsx=120, histnorm="probability density", name="samples", opacity=0.6),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(x=xgrid, y=loglaplace_pdf_std(xgrid, c_mc), mode="lines", name="theory"),
row=1,
col=1,
)
# Log-scale samples should look Laplace
fig.add_trace(
go.Histogram(x=z_samp, nbinsx=120, histnorm="probability density", name="log samples", opacity=0.6),
row=1,
col=2,
)
fig.add_trace(
go.Scatter(x=zgrid, y=laplace_pdf, mode="lines", name="Laplace pdf"),
row=1,
col=2,
)
fig.update_xaxes(title="x", type="log", row=1, col=1)
fig.update_xaxes(title="z = log x", row=1, col=2)
fig.update_yaxes(title="density", row=1, col=1)
fig.update_yaxes(title="density", row=1, col=2)
fig.update_layout(width=980, height=380)
fig.show()
# Monte Carlo check of mean/variance when they exist (c>2)
print("theory mean:", loglaplace_mean(c_mc))
print("MC mean :", x_samp.mean())
print("theory var :", loglaplace_variance(c_mc))
print("MC var :", x_samp.var())
theory mean: 1.1904761904761905
MC mean : 1.1881560637980781
theory var : 1.3605442176870748
MC var : 0.9558356968621015
9) SciPy Integration (scipy.stats.loglaplace)#
SciPy provides a ready-to-use distribution object:
pdf,logpdf,cdf,ppfrvsfor samplingfitfor MLE of parameters
We’ll use fit in two ways:
Shape-only fit with
loc=0, scale=1fixed (matches the “standard” model).Full
c, loc, scalefit (often less stable; interpret carefully).
from scipy.stats import loglaplace
c0 = 2.0
x = loglaplace.rvs(c0, size=5, random_state=rng)
print("SciPy rvs:", x)
# Evaluate
print("pdf:", loglaplace.pdf(x, c0))
print("cdf:", loglaplace.cdf(x, c0))
# Fit: generate data from the standard model
data = loglaplace_rvs_numpy(c0, size=5_000, rng=rng)
# Closed-form MLE for the standard model
c_hat_closed = loglaplace_mle_c(data)
# SciPy fit with loc/scale fixed to match the standard model
c_hat_scipy, loc_hat, scale_hat = stats.loglaplace.fit(data, floc=0.0, fscale=1.0)
print("c0 =", c0)
print("c_hat closed=", c_hat_closed)
print("c_hat SciPy =", c_hat_scipy)
print("(loc,scale) =", (loc_hat, scale_hat))
# Full fit (c, loc, scale) — can be sensitive for heavy tails
c_fit, loc_fit, scale_fit = stats.loglaplace.fit(data)
print("full fit (c,loc,scale)=", (c_fit, loc_fit, scale_fit))
SciPy rvs: [2.60896 0.72759 0.92112 0.64573 1.02715]
pdf: [0.05631 0.72759 0.92112 0.64573 0.92279]
cdf: [0.92654 0.26469 0.42424 0.20849 0.52608]
c0 = 2.0
c_hat closed= 2.004309521405735
c_hat SciPy = 2.004309521405735
(loc,scale) = (0.0, 1.0)
full fit (c,loc,scale)= (2.0137861206612757, -0.0034857193644913936, 1.000027126069282)
10) Statistical Use Cases#
10.1 Hypothesis testing#
In the standard model, inference on \(c\) is especially simple.
Since \(|\log X|\sim \mathrm{Exp}(\text{rate}=c)\), the statistic
satisfies
This gives exact confidence intervals and tests for \(c\).
10.2 Bayesian modeling#
The standard-model likelihood for \(c\) is proportional to
so a Gamma prior on \(c\) is conjugate:
(Here \(b\) is a rate parameter.)
10.3 Generative modeling#
If you want multiplicative heavy-tailed noise, you can model
which is equivalent to adding Laplace noise in log-space:
from scipy.stats import chi2
# Exact CI and an exact two-sided test for c in the standard model
c0 = 2.0
n = 300
x = loglaplace_rvs_numpy(c0, size=n, rng=rng)
S = np.sum(np.abs(np.log(x)))
alpha = 0.05
ci = (
chi2.ppf(alpha / 2, df=2 * n) / (2 * S),
chi2.ppf(1 - alpha / 2, df=2 * n) / (2 * S),
)
c_hat = loglaplace_mle_c(x)
test_stat = 2 * c0 * S # under H0, this is chi2_{2n}
p_left = chi2.cdf(test_stat, df=2 * n)
p_two_sided = 2 * min(p_left, 1 - p_left)
print("c0 =", c0)
print("c_hat=", c_hat)
print("95% exact CI for c:", ci)
print("two-sided exact p-value for H0: c=c0:", p_two_sided)
c0 = 2.0
c_hat= 2.0630561150287168
95% exact CI for c: (1.836183726795856, 2.3029522418961106)
two-sided exact p-value for H0: c=c0: 0.6061428639658745
# Conjugate Bayesian update for c (standard model)
# Prior: c ~ Gamma(a0, rate=b0)
a0, b0 = 2.0, 1.0
c0 = 2.0
n = 200
x = loglaplace_rvs_numpy(c0, size=n, rng=rng)
S = np.sum(np.abs(np.log(x)))
a_post = a0 + n
b_post = b0 + S
posterior_mean = a_post / b_post
posterior_ci = (
stats.gamma.ppf(0.025, a=a_post, scale=1 / b_post),
stats.gamma.ppf(0.975, a=a_post, scale=1 / b_post),
)
print("prior (a,b)=", (a0, b0))
print("posterior (a,b)=", (a_post, b_post))
print("posterior mean=", posterior_mean)
print("posterior 95% CI=", posterior_ci)
prior (a,b)= (2.0, 1.0)
posterior (a,b)= (202.0, 93.21649390394049)
posterior mean= 2.166998473555129
posterior 95% CI= (1.8784506146804234, 2.4758606604180544)
11) Pitfalls#
Invalid parameters: \(c\le 0\) (and
scale <= 0in SciPy) is invalid.Nonexistent moments: mean requires \(c>1\), variance requires \(c>2\), etc. For smaller \(c\), Monte Carlo estimates like the sample mean can look “unstable” because the target quantity is infinite.
Extreme samples: the inverse CDF produces very large values when \(U\) is extremely close to 1. In floating point code, clip uniforms away from 0/1.
Numerical evaluation: prefer
logpdffor likelihoods to avoid underflow, especially when \(x\) is tiny/huge or \(c\) is large.Fitting with
loc/scale: unconstrainedfitmay be sensitive for heavy-tailed data; consider fixingloc/scalewhen you have a clear standardization.
12) Summary#
Log-Laplace models positive, heavy-tailed data and is Laplace in log-space.
The PDF is piecewise power-law and the right tail behaves like a Pareto tail with index \(c\).
Raw moments satisfy \(\mathbb{E}[X^k]=\frac{c^2}{c^2-k^2}\) for \(|k|<c\), so some moments may not exist.
Sampling is easy with an inverse CDF (NumPy-only).
In the standard model, \(|\log X|\) is exponential, giving a closed-form MLE, exact chi-square inference, and a conjugate Gamma posterior for \(c\).